| Universidad | Universidad Nacional del Oeste |
| Carrera | Esp. en Ciencia de Datos |
| Materia | Fundamentos de Estadísticas (01050) |
| Profesora | Mg. Sc. Silvia N. Pérez |
| Grupo | Nº 3 |
| Alumnos | {Lic. Leticia Sosa, Ing. Federico Czerniawski, Mg. Pablo Pandolfo} |
| Fecha | Junio 2025 |
library(readxl)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(hrbrthemes)
library(plotly)
library(tidyr)
library(tidyverse)
library(scales)
library(corrplot)
library(car)
library(nortest)
library(tseries)
library(knitr)
library(vcd)
# Borramos ambiente de trabajo
rm(list = ls())
# Seteamos directorio de trabajo
setwd("/Users/ppando/Materias/data/materias/estadistica/tp_grupal")
# Cargamos datos del archivo
datos <- read_xlsx("seguros.xlsx")
# Mostramos los primeros 6 casos
head(datos)
# Mostramos estructura de los datos
str(datos)
## tibble [669 × 7] (S3: tbl_df/tbl/data.frame)
## $ Edad : num [1:669] 18 37 29 54 19 30 63 24 27 59 ...
## $ Genero : chr [1:669] "masculino" "femenino" "femenino" "femenino" ...
## $ IMC : num [1:669] 17.3 17.3 38.8 27.6 33.1 ...
## $ Fuma : chr [1:669] "si" "no" "no" "no" ...
## $ Region : chr [1:669] "nordeste" "nordeste" "sureste" "noroeste" ...
## $ ingreso: num [1:669] 1337 1181 1092 1348 1203 ...
## $ premio : num [1:669] 40.7 41 39.7 41.7 41.3 ...
# Contamos número de nulos por variables (columnas)
sapply(datos, function(x) sum(is.na(x)))
## Edad Genero IMC Fuma Region ingreso premio
## 0 0 0 1 0 0 0
# Eliminamos casos (filas) con nulos
datos <- datos[!is.na(datos$Fuma),]
# Verificamos que el caso con nulos fue eliminado
sapply(datos, function(x) sum(is.na(x)))
## Edad Genero IMC Fuma Region ingreso premio
## 0 0 0 0 0 0 0
seguros_long_Genero <- datos %>%
select(Edad, Genero, IMC, ingreso, premio) %>%
pivot_longer(cols = -Genero, names_to = "variable", values_to = "valor")
ggplotly(
ggplot(seguros_long_Genero, aes(x = Genero, y = valor, fill = Genero)) +
geom_boxplot(alpha = 0.7) +
facet_wrap(~ variable, scales = "free_y") + # Un panel por variable
theme_bw() +
theme(legend.position = "none") +
labs(
title = "Distribución de variables por Genero",
x = "Género",
y = "Valor"
)
)
seguros_long_Fuma <- datos %>%
select(Edad, IMC, Fuma, ingreso, premio) %>%
pivot_longer(cols = -Fuma, names_to = "variable", values_to = "valor")
ggplotly(
ggplot(seguros_long_Fuma, aes(x = Fuma, y = valor, fill = Fuma)) +
geom_boxplot(alpha = 0.7) +
facet_wrap(~ variable, scales = "free_y") + # Un panel por variable
theme_bw() +
theme(legend.position = "none") +
labs(
title = "Distribución de variables por Condición de Fumador",
x = "Fuma",
y = "Valor"
)
)
seguros_long_Region <- datos %>%
select(Edad, IMC, Region, ingreso, premio) %>%
pivot_longer(cols = -Region, names_to = "variable", values_to = "valor")
ggplotly(
ggplot(seguros_long_Region, aes(x = Region, y = valor, fill = Region)) +
geom_boxplot(alpha = 0.7) +
facet_wrap(~ variable, scales = "free_y") + # Un panel por variable
theme_bw() +
theme(legend.position = "none") +
labs(
title = "Distribución de variables por Region",
x = "Región",
y = "Valor"
)
)
Se observa que las variables cuantitativas agrupadas por las distintas cualitativas presentan distribuciones similares.
Además, la variable premio presenta un valor muy atípico.
Como se concluye del análisis anterior, la variable premio presenta un valor muy atípico.
Procedemos a eliminarlo.
gplot = datos %>% filter(Genero == "femenino") %>%
ggplot(mapping = aes(x=Region, y=premio, color=Fuma)) +
geom_point() +
stat_summary(fun ="mean", geom="crossbar", color="red", linewidth=0.5) + # mostramos la media
theme_bw() + xlab("Region") + ylab("Premio")
ggplotly(gplot)
# Lo eliminamos
datos <- datos %>% filter(premio != max(premio))
# Declaramos a las variables cualitativas como factor, para poder reconocer categorias
datos$Genero <- as.factor(datos$Genero)
datos$Fuma <- as.factor(datos$Fuma)
datos$Region <- as.factor(datos$Region)
# Mostramos medidas resumen
summary(datos)
## Edad Genero IMC Fuma Region ingreso
## Min. :18.00 femenino :321 Min. :15.96 no:535 nordeste:162 Min. : 417.7
## 1st Qu.:27.00 masculino:346 1st Qu.:26.18 si:132 noroeste:165 1st Qu.: 938.6
## Median :39.00 Median :30.21 sureste :168 Median :1087.3
## Mean :39.32 Mean :30.47 suroeste:172 Mean :1100.9
## 3rd Qu.:51.00 3rd Qu.:34.43 3rd Qu.:1263.8
## Max. :64.00 Max. :52.58 Max. :1756.6
## premio
## Min. :32.62
## 1st Qu.:37.48
## Median :38.84
## Mean :38.84
## 3rd Qu.:40.10
## Max. :44.67
# Volvemos a mostrar graficamente las medidas resumen de las variables cuantis agrupadas por las variables cualis
seguros_long_Genero <- datos %>%
select(Edad, Genero, IMC, ingreso, premio) %>%
pivot_longer(cols = -Genero, names_to = "variable", values_to = "valor")
ggplotly(
ggplot(seguros_long_Genero, aes(x = Genero, y = valor, fill = Genero)) +
geom_boxplot(alpha = 0.7) +
facet_wrap(~ variable, scales = "free_y") + # Un panel por variable
theme_bw() +
theme(legend.position = "none") +
labs(
title = "Distribución de variables por Genero",
x = "Género",
y = "Valor"
)
)
seguros_long_Fuma <- datos %>%
select(Edad, IMC, Fuma, ingreso, premio) %>%
pivot_longer(cols = -Fuma, names_to = "variable", values_to = "valor")
ggplotly(
ggplot(seguros_long_Fuma, aes(x = Fuma, y = valor, fill = Fuma)) +
geom_boxplot(alpha = 0.7) +
facet_wrap(~ variable, scales = "free_y") + # Un panel por variable
theme_bw() +
theme(legend.position = "none") +
labs(
title = "Distribución de variables por Condición de Fumador",
x = "Fuma",
y = "Valor"
)
)
seguros_long_Region <- datos %>%
select(Edad, IMC, Region, ingreso, premio) %>%
pivot_longer(cols = -Region, names_to = "variable", values_to = "valor")
ggplotly(
ggplot(seguros_long_Region, aes(x = Region, y = valor, fill = Region)) +
geom_boxplot(alpha = 0.7) +
facet_wrap(~ variable, scales = "free_y") + # Un panel por variable
theme_bw() +
theme(legend.position = "none") +
labs(
title = "Distribución de variables por Region",
x = "Región",
y = "Valor"
)
)
# Calculamos la matriz de correlación que nos muestra los coeficientes de correlación entre cada para de variables.
datos_numericos <- datos %>% select(Edad, IMC, ingreso, premio)
correlaciones <- cor(datos_numericos)
correlaciones
## Edad IMC ingreso premio
## Edad 1.000000000 0.14973714 0.008118905 0.01590152
## IMC 0.149737140 1.00000000 0.077763646 0.03828891
## ingreso 0.008118905 0.07776365 1.000000000 0.87995085
## premio 0.015901515 0.03828891 0.879950853 1.00000000
# Definimos una paleta de colores
colores <- colorRampPalette(c("blue", "white", "red"))(200)
# Visualizamos las correlaciones mediante un mapa de calor
corrplot(correlaciones, method = "color", col = colores, addCoef.col = "black")
# Visualizamos en un gráfico de dispersión la dependencia entre ingreso y premio
gplot = datos %>%
ggplot(mapping = aes(x=ingreso, y=premio, color=Genero)) +
geom_point() +
geom_smooth(formula= 'y ~ x', method = "lm", se = FALSE, color = "red") +
theme_bw() + xlab("Ingreso") + ylab("Premio") + ggtitle("Relación entre ingreso y premio")
ggplotly(gplot)
# Visualizamos en un gráfico de dispersión la dependencia entre Edad e IMC
gplot = datos %>%
ggplot(mapping = aes(x=Edad, y=IMC, color=Genero)) +
geom_point() +
geom_smooth(formula= 'y ~ x', method = "lm", se = FALSE, color = "red") +
theme_bw() + xlab("Edad") + ylab("IMC") + ggtitle("Relación entre Edad e IMC")
ggplotly(gplot)
En el gráfico se muestra claramente la correlación fuerte existente entre las variables ingreso y premio.
Además, las variables Edad e IMC presentan una correlación moderada.
Y las variables ingreso e IMC presentan una correlación débil.
El resto de las intersecciones no presentan correlaciones.
Ho -> No hay diferencias entre las medias del IMC de los grupos “Fuma = sí” y “Fuma = no”.
H1 -> Hay diferencias entre las medias de los dos grupos.
Tamaño de la muestra: 667
Nivel de Significancia: 0.05
Se consideran dos muestras independientes.
# Verificamos mediante un QQPlot
gplot <- datos %>%
ggplot(aes(sample = IMC, color = Fuma)) +
geom_qq() + geom_qq_line() +
facet_wrap(~ Fuma) +
ggtitle("QQ-plot de IMC por estado de fumador") +
theme_bw()
ggplotly(gplot)
Aparentemente ambas muestras siguen una distribución normal.
Verificamos con una función estadística (usamos Lilliefors porque el tamaño de la muestra es mayor a 50)
by(data = datos, INDICES = datos$Fuma, FUN = function(x) { lillie.test(x$IMC) })
## datos$Fuma: no
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$IMC
## D = 0.030947, p-value = 0.2444
##
## -------------------------------------------------------------------
## datos$Fuma: si
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$IMC
## D = 0.045829, p-value = 0.7123
p-value >= 0.05 en ambos grupos: NO hay evidencia estadística para rechazar H0, lo que indica que se puede suponer que los datos de la variable IMC siguen una distribución normal.
bartlett.test(IMC ~ Fuma, data=datos)
##
## Bartlett test of homogeneity of variances
##
## data: IMC by Fuma
## Bartlett's K-squared = 0.35468, df = 1, p-value = 0.5515
leveneTest(IMC ~ Fuma, data=datos) # Se usa cuando no se puede asegurar la normalidad
En ambos tests se observa un p-value alejado del 0.05, por lo cual no hay evidencia para rechazar H0, lo que permite suponer que hay homogeneidad de varianzas.
Como las muestras se suponen con distribución normal y las varianzas homogéneas, se utilizará el t.test para hacer la prueba de hipótesis.
t.test(IMC ~ Fuma, alternative = "two.sided", mu = 0, var.equal = TRUE, data=datos) # p-value = 0.2055
##
## Two Sample t-test
##
## data: IMC by Fuma
## t = -1.2673, df = 665, p-value = 0.2055
## alternative hypothesis: true difference in means between group no and group si is not equal to 0
## 95 percent confidence interval:
## -1.8990106 0.4092541
## sample estimates:
## mean in group no mean in group si
## 30.32296 31.06784
Como se pueden considerar ambas muestras con distribución normal, independientes y con varianzas homogéneas, se realizó un test t de Student para comparar las medias de IMC de los grupos “No fumadores” y “Fumadores”.
Como el p-valor está alejado y mayor del 0,05 NO hay evidencia suficiente para rechazar H0, lo que presupone que no existe una diferencia significativa en la variación del promedio del IMC entre ambos grupos de la población.
Ho -> No hay diferencias entre las medias del ingreso de los grupos “Genero = femenino” y “Genero = masculino”.
H1 -> Hay diferencias entre las medias de los dos grupos.
Tamaño de la muestra: 667
Nivel de Significancia: 0.05
Se consideran dos muestras independientes.
# Verificamos mediante un QQPlot
gplot <- datos %>%
ggplot(aes(sample = ingreso, color = Genero)) +
geom_qq() + geom_qq_line() +
facet_wrap(~ Genero) +
ggtitle("QQ-plot de ingreso por Género") +
theme_bw()
ggplotly(gplot)
Aparentemente ambas muestras siguen una distribución normal.
Verificamos con una función estadística (usamos Lilliefors porque el tamaño de la muestra es mayor a 50)
by(data = datos, INDICES = datos$Genero, FUN = function(x) { lillie.test(x$ingreso) })
## datos$Genero: femenino
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$ingreso
## D = 0.032072, p-value = 0.5839
##
## -------------------------------------------------------------------
## datos$Genero: masculino
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$ingreso
## D = 0.029988, p-value = 0.6313
p-value >= 0.05 en ambos grupos: NO hay evidencia estadística para rechazar H0, lo que indica que se puede suponer que los datos de la variable ingreso siguen una distribución normal cuando la agrupamos por Genero
bartlett.test(ingreso ~ Genero, data=datos)
##
## Bartlett test of homogeneity of variances
##
## data: ingreso by Genero
## Bartlett's K-squared = 6.4844, df = 1, p-value = 0.01088
leveneTest(ingreso ~ Genero, data=datos) # Se usa cuando no se puede asegurar la normalidad
En ambos tests se observa un p-value menor al 0.05, por lo cual se rechaza H0, que se interpreta que las varianzas no son homogéneas.
Como las muestras se suponen con distribución normal y las varianzas heterogénas, se utilizará el test de Welch para hacer la prueba de hipótesis.
oneway.test(ingreso ~ Genero, data=datos)
##
## One-way analysis of means (not assuming equal variances)
##
## data: ingreso and Genero
## F = 161.87, num df = 1.00, denom df = 662.21, p-value < 2.2e-16
t.test(ingreso ~ Genero, alternative = "two.sided", mu = 0, var.equal = FALSE, data=datos)
##
## Welch Two Sample t-test
##
## data: ingreso by Genero
## t = -12.723, df = 662.21, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group femenino and group masculino is not equal to 0
## 95 percent confidence interval:
## -245.7775 -180.0568
## sample estimates:
## mean in group femenino mean in group masculino
## 990.4566 1203.3738
En ambos test se obtuvo un p-valor muy pequeño lo que indica que se rechaza H0, por lo tanto se puede decir que, estadísticamente, la media poblacional de ingresos entre ambos géneros es diferente.
t.test(ingreso ~ Genero, alternative = "greater", mu = 0, var.equal = FALSE, data=datos)
##
## Welch Two Sample t-test
##
## data: ingreso by Genero
## t = -12.723, df = 662.21, p-value = 1
## alternative hypothesis: true difference in means between group femenino and group masculino is greater than 0
## 95 percent confidence interval:
## -240.4826 Inf
## sample estimates:
## mean in group femenino mean in group masculino
## 990.4566 1203.3738
Conclusión: Se puede observar que el p-valor = 1 por lo que se NO se rechaza la Hipótesis Nula, confirmando estadísticamente que la media poblacional de ingresos del género masculino es mayor a la media de ingresos del genero femenino.
Ho -> No hay diferencias entre las medias de IMC, ingreso y premio por Region.
H1 -> Hay alguna media diferente entre las regiones.
Tamaño de la muestra: 667
Nivel de Significancia: 0.05
Se consideran cuatro muestras independientes.
Ahora vamos a analizar IMC por Región
# Verificamos mediante un QQPlot
gplot <- datos %>%
ggplot(aes(sample = IMC, color = Region)) +
geom_qq() + geom_qq_line() +
facet_wrap(~ Region) +
ggtitle("QQ-plot de IMC por Región") +
theme_bw()
ggplotly(gplot)
Aparentemente las muestras siguen una distribución normal.
Verificamos con una función estadística (usamos Lilliefors porque el tamaño de la muestra es mayor a 50)
by(data = datos, INDICES = datos$Region, FUN = function(x) { lillie.test(x$IMC) })
## datos$Region: nordeste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$IMC
## D = 0.04531, p-value = 0.5743
##
## -------------------------------------------------------------------
## datos$Region: noroeste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$IMC
## D = 0.06141, p-value = 0.134
##
## -------------------------------------------------------------------
## datos$Region: sureste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$IMC
## D = 0.040232, p-value = 0.7278
##
## -------------------------------------------------------------------
## datos$Region: suroeste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$IMC
## D = 0.040478, p-value = 0.7025
p-value >= 0.05 en los cuatro grupos: NO hay evidencia estadística para rechazar H0, lo que indica que se puede suponer que los datos de la variable IMC siguen una distribución normal cuando la agrupamos por Región
bartlett.test(IMC ~ Region, data=datos)
##
## Bartlett test of homogeneity of variances
##
## data: IMC by Region
## Bartlett's K-squared = 9.978, df = 3, p-value = 0.01875
leveneTest(IMC ~ Region, data=datos) # Se usa cuando no se puede asegurar la normalidad
En ambos tests se observa un p-value menor al 0.05, por lo cual se rechaza H0, por lo que se interpreta que las varianzas no son homogéneas.
Como las muestras se suponen con distribución normal y las varianzas heterogénas, se utilizará el test de Welch para hacer la prueba de hipótesis.
oneway.test(IMC ~ Region, data=datos)
##
## One-way analysis of means (not assuming equal variances)
##
## data: IMC and Region
## F = 19.231, num df = 3.00, denom df = 367.54, p-value = 1.325e-11
kruskal.test(IMC ~ Region, data = datos)
##
## Kruskal-Wallis rank sum test
##
## data: IMC by Region
## Kruskal-Wallis chi-squared = 51.743, df = 3, p-value = 3.397e-11
En ambos tests se obtuvo un p-valor muy pequeño lo que indica que se rechaza H0, por lo tanto se puede decir que, estadísticamente, la media poblacional de IMC entre las regiones es diferente.
Se va a utilizar una función para evaluar de a pares de regiones las medias de IMC, a fin de identificar cuales son las medias que difieren.
comparar_dos_regiones = function(ds, nombre_variable){
regiones_unicas = ds %>% distinct(Region)
lst = as.list(as.data.frame(t(regiones_unicas)))
pvalores_ds = data.frame(variables = character(), p_valor=character())
formula_str <- reformulate("Region", nombre_variable)
while(length(lst)>1){
region_A = as.character(lst[1])
lst = lst[-1]
for(region_B in lst){
region_B = as.character(region_B)
working_ds = ds %>% filter(Region %in% c(region_A, region_B))
rdo = t.test(formula_str, alternative = "two.sided", mu = 0, var.equal = FALSE, data=working_ds)
label = paste0(region_A,"-",region_B)
pvalores_ds = pvalores_ds %>% union(data.frame(variables=label, p_valor=as.character(rdo[["p.value"]])))
}
}
return(pvalores_ds)
}
dt = comparar_dos_regiones(datos, "IMC")
kable(dt, caption="<center><b>Resultado Test t de Student (p-valor) - IMC cada 2 Regiones</b></center>")
Table:
| variables | p_valor |
|---|---|
| nordeste-sureste | 2.29543858727147e-11 |
| nordeste-noroeste | 0.537331304440342 |
| nordeste-suroeste | 0.00171066577169515 |
| sureste-noroeste | 4.13625651114826e-10 |
| sureste-suroeste | 7.79698255660573e-05 |
| noroeste-suroeste | 0.00993143266327661 |
De acuerdo al test t de Student, todos tienen los p-valores muy cercanos a 0, por lo que se rechaza H0 en todos los casos excepto para el par de regiones “nordeste-noroeste”.
Analizando el par “nordeste-noroeste” se puede observar que el p-valor es 0.5373, muy alto, por lo que no hay evidencia para rechazar H0 y se podria contemplar que estadísticamente la media de IMC de ambas poblaciones son iguales.
Ahora vamos a analizar ingreso por Región
# Verificamos mediante un QQPlot
gplot <- datos %>%
ggplot(aes(sample = ingreso, color = Region)) +
geom_qq() + geom_qq_line() +
facet_wrap(~ Region) +
ggtitle("QQ-plot de ingreso por Región") +
theme_bw()
ggplotly(gplot)
Aparentemente las muestras siguen una distribución normal.
Verificamos con una función estadística (usamos Lilliefors porque el tamaño de la muestra es mayor a 50)
by(data = datos, INDICES = datos$Region, FUN = function(x) { lillie.test(x$ingreso) })
## datos$Region: nordeste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$ingreso
## D = 0.048198, p-value = 0.4735
##
## -------------------------------------------------------------------
## datos$Region: noroeste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$ingreso
## D = 0.043916, p-value = 0.6098
##
## -------------------------------------------------------------------
## datos$Region: sureste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$ingreso
## D = 0.058215, p-value = 0.1775
##
## -------------------------------------------------------------------
## datos$Region: suroeste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$ingreso
## D = 0.043712, p-value = 0.5843
p-value >= 0.05 en los cuatro grupos: NO hay evidencia estadística para rechazar H0, lo que indica que se puede suponer que los datos de la variable ingreso siguen una distribución normal cuando la agrupamos por Región
bartlett.test(ingreso ~ Region, data=datos)
##
## Bartlett test of homogeneity of variances
##
## data: ingreso by Region
## Bartlett's K-squared = 8.4901, df = 3, p-value = 0.0369
leveneTest(ingreso ~ Region, data=datos) # Se usa cuando no se puede asegurar la normalidad
En ambos tests se observa un p-value menor al 0.05, por lo cual se rechaza H0, por lo que se interpreta que las varianzas no son homogéneas.
Como las muestras se suponen con distribución normal y las varianzas heterogénas, se utilizará el test de Welch para hacer la prueba de hipótesis.
oneway.test(ingreso ~ Region, data=datos)
##
## One-way analysis of means (not assuming equal variances)
##
## data: ingreso and Region
## F = 1.164, num df = 3.00, denom df = 367.84, p-value = 0.3233
kruskal.test(ingreso ~ Region, data = datos)
##
## Kruskal-Wallis rank sum test
##
## data: ingreso by Region
## Kruskal-Wallis chi-squared = 2.9365, df = 3, p-value = 0.4015
En ambos tests se obtuvo un p-valor > 0.05 lo que indica que NO hay evidencia para rechazar H0, por lo tanto se puede decir que, estadísticamente, la media poblacional de ingreso entre las regiones es similar.
Se va a utilizar una función para evaluar de a pares de regiones las medias de ingreso, a fin de corroborar la conclusión anterior.
dt = comparar_dos_regiones(datos, "ingreso")
kable(dt, caption="<center><b>Resultado Test t de Student (p-valor) - ingreso cada 2 Regiones</b></center>")
Table:
| variables | p_valor |
|---|---|
| nordeste-sureste | 0.228525487375051 |
| nordeste-noroeste | 0.608749809791648 |
| nordeste-suroeste | 0.526807312402343 |
| sureste-noroeste | 0.0792080745348705 |
| sureste-suroeste | 0.604634481804345 |
| noroeste-suroeste | 0.255010367991364 |
Ahora vamos a analizar premio por Región
# Verificamos mediante un QQPlot
gplot <- datos %>%
ggplot(aes(sample = premio, color = Region)) +
geom_qq() + geom_qq_line() +
facet_wrap(~ Region) +
ggtitle("QQ-plot de Premio por Región") +
theme_bw()
ggplotly(gplot)
Aparentemente las muestras siguen una distribución normal.
Verificamos con una función estadística (usamos Lilliefors porque el tamaño de la muestra es mayor a 50)
by(data = datos, INDICES = datos$Region, FUN = function(x) { lillie.test(x$premio) })
## datos$Region: nordeste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$premio
## D = 0.045881, p-value = 0.554
##
## -------------------------------------------------------------------
## datos$Region: noroeste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$premio
## D = 0.046633, p-value = 0.5126
##
## -------------------------------------------------------------------
## datos$Region: sureste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$premio
## D = 0.052678, p-value = 0.3056
##
## -------------------------------------------------------------------
## datos$Region: suroeste
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$premio
## D = 0.029289, p-value = 0.9742
p-value >= 0.05 en los cuatro grupos: NO hay evidencia estadística para rechazar H0, lo que indica que se puede suponer que los datos de la variable premio siguen una distribución normal cuando la agrupamos por Región
bartlett.test(premio ~ Region, data=datos)
##
## Bartlett test of homogeneity of variances
##
## data: premio by Region
## Bartlett's K-squared = 3.4976, df = 3, p-value = 0.3211
leveneTest(premio ~ Region, data=datos) # Se usa cuando no se puede asegurar la normalidad
En ambos tests se observa un p-value mayor a 0.05, por lo cual NO hay evidencia para rechazar H0, por lo que se interpreta que las varianzas son homogéneas.
Como las muestras se suponen con distribución normal, las varianzas son homogéneas y se trata de varios grupos, se utilizará el test ANOVA para hacer la prueba de hipótesis.
anova_result <- aov(premio ~ Region, data = datos)
summary(anova_result)
## Df Sum Sq Mean Sq F value Pr(>F)
## Region 3 2.4 0.801 0.187 0.905
## Residuals 663 2833.3 4.273
oneway.test(premio ~ Region, data=datos, var.equal = TRUE) # p-value = 0.905
##
## One-way analysis of means
##
## data: premio and Region
## F = 0.18744, num df = 3, denom df = 663, p-value = 0.905
En ambos tests se obtuvo un p-valor muy alto lo que indica que hay una fuerte evidencia para NO rechazar H0, por lo tanto se puede decir que, estadísticamente, la media poblacional de premio entre las regiones es similar.
Se va a utilizar una función para evaluar de a pares de regiones las medias de premio, a fin de corroborar la conclusión anterior.
dt = comparar_dos_regiones(datos, "premio")
kable(dt, caption="<center><b>Resultado Test t de Student (p-valor) - Premio cada 2 Regiones</b></center>")
Table:
| variables | p_valor |
|---|---|
| nordeste-sureste | 0.688064368252902 |
| nordeste-noroeste | 0.732396906064601 |
| nordeste-suroeste | 0.925471665312624 |
| sureste-noroeste | 0.445882813507527 |
| sureste-suroeste | 0.627025171971952 |
| noroeste-suroeste | 0.813280564328186 |
Procedemos a ejecutar dos tests Post-HOC: test de Tukey y Bonferroni
intervalos = TukeyHSD(anova_result)
plot(intervalos)
pairwise.t.test(x = datos$premio, g = datos$Region, p.adjust.method = "bonferroni",
pool.sd = TRUE, paired = FALSE, alternative = "two.sided")
##
## Pairwise comparisons using t tests with pooled SD
##
## data: datos$premio and datos$Region
##
## nordeste noroeste sureste
## noroeste 1 - -
## sureste 1 1 -
## suroeste 1 1 1
##
## P value adjustment method: bonferroni
Con ambos tests corroboramos la conclusión del ANOVA
Para evaluar dependencia vamos a realizar Tablas de Contingencia, que nos permiten evaluar dos variables cualitativas y buscar una relación entre ellas.
Comparamos la tabla de frecuencias observadas con la tabla de frecuencias esperada (suponiendo independencia) utilizando un estadístico de prueba.
Siendo >> X = Región & Y = Fuma
frecuencias = table(datos$Region, datos$Fuma)
margin.table(frecuencias, 1)
##
## nordeste noroeste sureste suroeste
## 162 165 168 172
margin.table(frecuencias, 2)
##
## no si
## 535 132
prop.table(frecuencias)
##
## no si
## nordeste 0.19640180 0.04647676
## noroeste 0.19790105 0.04947526
## sureste 0.19040480 0.06146927
## suroeste 0.21739130 0.04047976
chisq.test(frecuencias, correct = F)
##
## Pearson's Chi-squared test
##
## data: frecuencias
## X-squared = 4.1168, df = 3, p-value = 0.2491
Conclusión: En el test de Chi cuadrado sobre la tabla de frencuencias se observa un p-valor 0.2491, es bastante mayor a 0.05, no hay evidencia para rechazar H0.
Esto indicaría que estadísticamente se puede considerar que hay independencia entre ambas variables cualitativas: Region - Fuma
Entonces, NO hay evidencia suficiente para afirmar que la condición de fumador dependa de la región.
assocstats(frecuencias)
## X^2 df P(> X^2)
## Likelihood Ratio 4.1061 3 0.25023
## Pearson 4.1168 3 0.24913
##
## Phi-Coefficient : NA
## Contingency Coeff.: 0.078
## Cramer's V : 0.079
Corroborando los niveles de asociación se obtienen indices muy bajos lo que refuerza la evidencia de independencia obtenida por el test de Chi cuadrado.
FIN